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ABSTRACT 

We critically examine the reliability of the dark matter clustering in high-resolution cosmological N- 
body simulations on scales below the mean separation length of particles. The particle discreteness effect 
imposes the two fundamental limitations on those scales; the lack of the initial fluctuation power and 
the finite mass resolution. We address this problem applying the dark halo approach and are able to 
discuss separately how those two limitations affect the dark matter clustering in Y-body simulations at 
early epochs. We find that limitations of the dark matter clustering are primarily determined by the 
mass of particles. By a detailed comparison with three major cosmological simulations, we also find that 
in order to reproduce a proper amplitude of the dark matter clustering on small scales, halos with a 
characteristic nonlinear mass, M-^<ii^{z) defined by gr^M^i^^z) = I, must be resolved in the simulation. 
This leads to a critical redshift Zcrit determined by MNL(^crit) = "-haio'Tipart where rihaio is the number 
of particles necessary to resolve the typical nonlinear mass halo (~ 10). We conclude that, at least as 
far as the two-point correlation functions are concerned, the dark matter clustering in high-resolution 
Y-body simulations on scales below the mean particle separation is reliable down to the gravitational 
force resolution length only for z < Zcrit, while it is strongly affected by the finite mass resolution for 

^ ^ ^crit • 

Subject headings: gravitation — cosmology: theory - dark matter ~ large-scale structure of universe - 
methods: Y-body simulations 



1. INTRODUCTION 

Y-body simulations have played a key role in study- 
ing large-scale structure formation in the universe. In 
particular, they enable one to probe the quasi-linear to 
nonlinear clustering in a direct manner without resorting 
to any unrealistic approximations. Thus they are exten- 
sively used in cosmology as a powerful and indispensable 
tool complementary to analytical techniques. Cosmologi- 
cal Y-body simulations have been first applied in explor- 
ing the nonlinear growth of two-point correlation functions 
(Miyoshi & Kihara 1975; Aarseth, Gott, & Turner 1979), 
and in studying the halo mass function (Press & Schechter 
1974). A seminal paper by Davis et al. (1985) clearly 
demonstrated that this methodology provides a variety of 
testable predictions for a given cosmological model. Since 
then almost all precise predictions in gravitational non- 
linear phenomena have relied heavily on N-body simula- 
tions, including accurate nonlinear modeling of the dark 
matter power spectrum (Hamilton et al. 1991; Peacock & 
Dodds 1996), redshift-space distortion (Suto & Suginohara 
1991; Cole, Fisher & Weinberg 1994; Magira, Jing & Suto 
2000), three- and four-point correlation functions (Suto 
1993; Matsubara & Suto 1994; Suto & Matsubara 1994), 
probability distribution function of the nonlinear density 
field (Suto, Itoh & Inagaki 1990; Lahav et al. 1993; Kayo, 



Taruya & Suto 2001), the mass function of virialized dark 
halos (Bahcall & Gen 1993; Sheth & Tormcn 1999; Jenk- 
ins et al. 2001), the density profiles of dark matter halos 
(Navarro, Frenk & White 1996,1997; Fukushige & Makino 
1997; Moore et al. 1998; Jing & Suto 2000; Yoshida et al. 
2000), dark halo biasing (Watanabe, Matsubara, & Suto 
1994; Jing & Suto 1998; Jing 1998; Sheth & Tormen 1999; 
Colberg et al. 2000; Taruya et al. 2001), and clustering 
on the past light-cone (Hamana et al. 2001a, b; Evrard et 
al. 2001). 

Of course, the reliability of Y-body simulations is 
severely limited by various numerical restrictions. Among 
them, the effective dynamical range is quite important. 
Consider a simulation employing Ypart particles in a co- 
moving cubic box of i^ox- Obviously the mass resolution 
is set by the number of particles, and the corresponding 
resolution in length scale is up to the mean particle sep- 

1/3 

aration Apart — ^box/Ypart- In fact, the initial perturba- 
tions of simulations are usually set on uniform grid points 
of an interval Apart- Therefore one cannot properly set 
the desired initial power spectrum for k > fcNyquist where 
^Nyquist = '''/Apart IS the Nyquist wavenumber^. Neverthe- 
less most N-body codes adopt a force resolution (e.g., the 
gravitational softening length egrav) typically one order of 
magnitude smaller than Apart- Are these high-resolution 



This is also the case for a Af-body simulation that adopts a "glass" distribution for particles to represent an initial uniform state (White 
1996, Jenkins et al. 1998), because the Nyquist wavenumber in this case is also set by ^Nyquist = ""/Apart. 

^ Throughout the present paper, we use the term "high-resolution simulation" to refer to the iV-body codes that employ the gravitational 
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simulations ^ really reliable on scales below the particle 
mean separation length Apart where the initial conditions 
cannot be set properly ? Although several attempts have 
been made to this direction (e.g., Efstathiou et al. 1985; 
1988; Gclb & Bcrtschingcr 1994; Mclott ct al. 1997; Splin- 
ter et al. 1998), this natural question has remained unan- 
swered in a quantitative manner, and is exactly what we 
address in detail in the rest of the paper. 

Qualitatively speaking, high-resolution simulations are 
justified only if almost all the power on scales below the 
mean separation Apart later epochs is generated via non- 
linear mode- coupling from the power on larger scales. The 
transfer of the fluctuation power in cosmology due to the 
nonlinear mode-coupling is examined numerically by Sugi- 
nohara ct al. (1991), as well as by perturbation analysis in 
Suto & Sasaki (1991) and Makino, Sasaki & Suto (1992). 
According to their results, the nonlinear power transfer is 
expected to eventually dominate the initial power if the ef- 
fective spectral index at the scale, ricH' ^ rflogP(fc)/(ilogfc, 
is much less than —1. Actually this is indeed the case for 
most realistic cosmological scenarios including cold dark 
matter (CDM) models whose power spectrum asymptoti- 
cally approaches oc k^^. Nevertheless the above qualita- 
tive argument does not answer precisely when the results 
below Apart become reliable. 

To be more specific, we examine the evolution of dark 
matter clustering in high-resolution cosmological A''-body 
simulations taking explicit account of the two fundamental 
limitations: (i) the lack of the initial fluctuation power on 
scales smaller than the mean particle separation, and (ii) 
the finite mass resolution. To do this, we apply the dark 
halo approach in which the nonlinear dark matter two- 
point correlation function is described by a sum of two 
contributions; one from the correlations between the ha- 
los (halo-halo correlation) , and the other from correlations 
of matter within the same halo (the Poisson contribution) 
(Peebles 1980 and references therein, for recent develop- 
ments in the dark halo approach, see Ma & Fry 2000; Sel- 
jak 2000). We show in the following sections that this 
approach provides a simple but powerful way to examine 
separately how those two limitations affect the simulation 
results. This approach is quite different from previous ones 
which usually attempted to find convergence among dif- 
ferent N-body codes with different numerical resolutions 
(e.g., Efstathiou et al. 1985, 1988; Gelb & Bertschinger 
1994; Melott ct al. 1997; SpUntcr ct al. 1998; Moore et al. 
1998; Jenkins et al. 1998, 2001; Knebe et al. 2000). 

We focus on the two-point correlation function of dark 
matter to examine the limitations in high-resolution A'^- 
body simulations for the following reasons. Firstly, the 
two-point correlation function is one of the most funda- 
mental statistics to quantify clustering properties, and has 
been studied in detail both analytically and numerically. 
Actually, there exist accurate fitting formulae of the dark 
matter power spectrum (Hamilton et al. 1991; Peacock 
& Dodds 1996) that reproduce the behavior of the two- 
point correlation function in linear to highly nonlinear 
regimes. Secondly, the estimator of the two-point corre- 
lation function from the particle distribution in TV-body 
data is well established (Kerscher, Szapudi & Szalay 2000 
and references therein). These two advantages enable us 

softening length smaller than the mean particle separation such as I 



to investigate limitations in high-resolution A'^-body sim- 
ulations through detailed comparison between numerical 
results and model predictions. 

The outline of this paper is as follows. In §2, we summa- 
rize our dark halo approach, and demonstrate the extent to 
which the two limitations may affect the evolution of dark 
matter clustering on scales below the mean particle sepa- 
ration length. We present detailed comparison between a 
high-resolution A''-body data and our model predictions in 
§3, and show that our model reproduces well an apparent 
depression in the dark matter two-point correlation func- 
tion measured from those simulation data. In particular 
we find that the depression becomes less appreciable as the 
ratio of the characteristic nonlinear mass Mnl to the parti- 
cle mass JTipart increases. On the basis of this, we evaluate 
the critical redshift ^crit beyond which the correlation func- 
tion of the dark matter in high-resolution A'-body simula- 
tions on scales below the mean particle separation length 
suffers significantly from discreteness limitations. Finally 
§4 is devoted to a summary and further discussion. 

2. DARK HALO APPROACH 

The dark halo approach has a long history (Neyman & 
Scott 1952; Limber 1953; Peebles 1974, 1980; McClelland 
& Silk 1977), and has recently applied to various problems 
in galaxy formation and cosmological nonlinear cluster- 
ing (Scherrer & Bertschinger 1991; White & Frenk 1991; 
Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; 
Sheth & Jain 1997; Komatsu & Kitayama 1999; Somerville 
& Primack 1999; Seljak 2000; Peacock & Smith 2000; Ma 
& Pry 2000; Cooray, Hu & Miralda-Escude 2000; Cooray 
& Hu 2001). In this section, we summarize several ex- 
pressions which are most relevant to the current analysis. 
Further details of this approach can be found in those pa- 
pers. 

2.1. Ingredients in m,odel predictions 

The accurate description of the matter clustering on 
the basis of the dark halo approach requires three major 
ingredients; the mass function of dark halos, the mass- 
dependent biasing factor of halos, and the density profile 
of halos. We present a summary of those three ingredients 
below. 

We adopt an analytical model of the halo mass function 
proposed by Sheth & Tormen (1999): 

±^dM^^fii.)di^ 
dM M'^^ ' 

= gA[l + (a.)-]x^^exp(-f)^, (1) 

where 



po is the mean cosmic mass density, and the numerical 
coefficients a and p are empirically fitted from N-body 
simulations as o = 0.707 and p = 0.3. Here cr(M, z) is 
the root-mean-square fluctuations in the matter density 
top-hat smoothed over a scale Rm = (3M/47r/9o)^^'^ and 
5c{z) is the threshold over-density for spherical collapse 

and Tree codes. 



3 



(see Nakamura & Suto 1997 and Henry 2000 for useful fit- 
ting functions) . The normalization constant A in equation 
(1) is determined by requiring that 



J_ f —MdM = f fMdu = 1. 
po J dM J ' 



(3) 



A simple analytic model of the halo bias was proposed 
by Mo & White (1996) applying the extended Press- 
Schechter theory (e.g., Lacey & Cole 1993). Jing (1998) 
and Sheth & Tormen (1999) discussed the correction for 
the mass-dcpendcncc, and the stochastic nature of the 
halo biasing is considered subsequently by Taruya & Suto 
(2000) and Somerville et al. (2001). For the present pur- 
pose, we use the scale-independent model by Sheth & Tor- 
men (1999) which reasonably reproduces iV-body data: 



b{u) = 1 + 



au 



1 



+ 



2p 



5c{l + av) 



(4) 



Finally we have to specify the density profile of dark mat- 
ter halos. We adopt the NFW profile (Navarro, Frenk & 
White 1996): 



p{r) 



(r/r«)(l+r/r,)2' 



(5) 



where pg and denote the scaling density and radius. 
While more recent simulations (Fukushige & Makino 1997; 
Moore et al. 1998) suggest that the inner profile is indeed 
steeper, we adopt the original profile just for simplicity. It 
is conventional to introduce the concentration parameter 
c = rvir/^s, where Tvir is the virial radius of the halo. In 
terms of this, the mass of halo is written as 



M = 



iwpsr: 



,3 

vir 



log(l + c) - 



1 



(6) 



Since the spherical collapse model indicates that M = 
47rr^jj.(5vir(-2)/Oo/3 with 6vir{z) being the over-density of col- 
lapse (see, Nakamura & Suto 1997; Henry 2000), one can 
express ps in terms of M and c using equation (6) . Thus 
our description of dark halos is completed if the concen- 
tration parameter c is expressed as a function of M and z. 
We adopt the following empirical relations: 



c(M, z) 



l + z 

5 



M 



1O14/1-1M0 

M 



l + z\ 10"ft-iMp 



-0.13 



-0.13 



for ACDM, 



for SCDM. 



(7) 

The above c-M relation for ACDM model is suggested by 
Bullock et al. (2001), while for SCDM model, we adopt the 
smaller normalization found by Bartelmann et al. (1998). 
For simplicity, we do not take into account the scatter in 
c-M relation (Bullock et al. 2001), but it should have only 
secondary importance on a statistical measurement such 
like the two-point correlation function of the dark mat- 
ter. We compare the dark matter power spectrum pre- 
dicted by the dark halo approach (see next subsection) 
with that computed by the nonlinear fitting formula of 
Peacock & Dodds (1996). We find a good agreement be- 
tween them for the SCDM and ACDM models over the 
range 0.01 < fc[/iMpc~^] < 100 and over < z < 5. 
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Fig. 1. — Comparison between dark matter two-point correlation 
functions predicted from the dark halo approach with and without 
taking into account the two limitations in the A'^-body simulation. 
Left panels: The effect of the lack of the initial perturbation power 
on small scales is demonstrated. The dotted line shows the nonlin- 
ear correlation function computed using Peacock & Dodds (1996) 
fitting function. The solid line shows the two-point correlation func- 
tion from the halo-halo contribution without taking into account the 
small-scale cutoff, while the dashed lines from upper to lower are for 
Rmin = 27r/fcmax = 0.3, 1, 3/i~^Mpc. Top and bottom panels are 
for 2 = and z = 3, respectively. Right panels: The effect of the 
finite mass resolution is demonstrated. The solid line shows the 
correlation function from the Poisson contribution in the case of no 
low mass cutoff, while the dashed lines from upper to lower are for 
Mmin = 1, 3, 10x10"/i-1Mq. 



2.2. Predicting two-point correlation functions 

In the dark halo approach, the two-point correlation 
function of the dark matter is the sum of two contribu- 
tions; ^(r) = ^^^{r) -\- ^P(r). The first term ("the halo- 
halo contribution") arises due to the correlation between 
two halos of density profile p{r) that are biased tracers of 
the dark matter distribution. In Fourier space, this term 
can be explicitly written as (Seljak 2000): 

2 



P'^^ik) = Pii„(fc) 



J dyf{y)h{y)y{k;M) 



(8) 



where P\m{k) is the linear power spectrum, y(fc; M) is the 
Fourier transform of the halo profile normalized by its 
mass, y{k; M) = p{k,M)/M, and M is related to ly via 
equation (2). The second "Poisson term" is due to the 
correlation between dark matter particles in the same halo 
and is expressed as 

M{iy) 



P^(fc)= di^fiiy) 



Po 



\yik;M)\' 



(9) 



Then the dark matter two-point correlation function is 
simply given by the sum of two terms 

where X stands for either hh or P for the halo-halo and 
Poisson terms, respectively. 



1 dkk^ 


sin(fcr) 




kr 



P^ik). (10) 
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by the halo-halo contribution only, while the small-scale 
power is dominated by the Poisson contribution. 

To simulate the effect of the numerical limitations more 
quantitatively, we artificially set 

P^^{k) = for k>k^^^ (11) 



for the halo-halo term, and 



dn 
dM 



for M < M„ 



(12) 



0.1 1 

r [h"^Mpc] 

Fig. 2. — The two-point correlation function measured from the 
Virgo SCDM Af-body simulation (filled circles) are compared with 
predictions. Dashed lines shows the nonlinear prediction by Peacock 
& Dodds (1996). Predictions by the dark halo approach, to which 
corrections for the lack of small-scale perturbation power and the 
finite particle mass are made (see §2), are shown by dotted, long- 
dashed and solid lines for the halo-halo term, Poisson term and 
their sum, respectively. The arrows indicate the mean separation 
lengths Apart of simulation particles. The redshift of the realization 
is indicated in each panel. 

2.3. Effects of the lack of small-scale power and the finite 
mass resolution 

Before presenting the detailed comparison between sim- 
ulations and our model predictions, it is instructive to 
show the main features that we expect in the present 
methodology. Figure 1 plots the dark matter two-point 
correlation functions at ^ = ( Upper panels) and z = 3 
{Lower panels) for the set of parameters employed for the 
Virgo ACDM simulation (sec Tabic 1). Solid lines in the 
left and right panels show the halo-halo and the Poisson 
terms, respectively. For comparison, dotted lines indicate 
the full nonlinear two-point correlation functions of the 
dark matter computed with the Peacock & Dodds (1996) 
formula. Clearly the large-scale power is well described 

^ For details of the output formats of the Hubble volume simulation, see the Hubble volume project home page: 
garching.mpg.de/Virgo/hubble.html 



for the Poisson term. The latter is equivalent to integrat- 
ing equation (9) only for v > j/min = [Sc{z)/a{M^in, z)]'^. 
The condition (11) corresponds to the lack of the ini- 
tial fluctuation power beyond the Nyquist wavenumber 
^Nyquist ~ ^max, whilc the coudition (12) reflects the flnite 
mass resolution corresponding to the simulation particles. 
In principle, the condition (12) also affects the halo-halo 
term, equation (8), however we do not include it, be- 
cause our aim here is to extract the most essential feature 
considering the most important defect in the numerical 
simulation for each term. We have found, however, that 
including the correction results in only a minor change 
in the estimates for the P'*'' term, so we do not consider 
the condition (12) in the halo-halo term below. Specifi- 
cally, the dashed lines plotted in left and right panels of 
Figure 1 adopt i?,nin = T^/kmax ~ 0.3, 1, 3/i^^Mpc, and 
Mmin = lO"ft-iM0, 3 X lO"/i-iM0, IO^^H-^Mq, from 
top to bottom. 

The message from these plots is that the dark matter 
clustering below the mean separation length is really sen- 
sitive to the two parameters fcmax and Mi„i„. In particular, 
the effect is quite substantial at higher redshifts where the 
major fraction of mass is contained in less massive ha- 
los that cannot be resolved with a given particle mass in 
simulations. In the next section we compare our model 
predictions with the evolution of simulation data. 

3. COMPARISON WITH n-BODY DATA 

3.1. N-body data 

For the present purpose, we need the current best 
simulation data. So we decided to use VIRGO-SCDM 
and VIRGO-ACDM simulations (Jenkins et al. 1998) and 
HUBBLE-ACDM (Jenkins et al. 2001; Evrard et al. 2001), 
all of which are publicly available. The simulation param- 
eters are summarized in Table 1. Their initial conditions 
were generated using the CDM transfer functions com- 
puted by Bond & Efstathiou (1984) for the VIRGO sim- 
ulations, while CMBFAST (Seljak & Zaldarriaga 1996) is 
used for the Hubble volume simulation which is well ap- 
proximated by the analytic fitting function by Eisenstein & 
Hu (1998). For the VIRGO simulations, we use snapshot 
data at 2 = 1, 2 and 3 for SCDM, and at ^ = 1, 3 and 5 for 
ACDM. In the case of the Hubble volume ACDM simula- 
tion, we do not use the full snapshot data at z, but rather 
"deep wedge" light-cone output. This dataset subtends an 
81.45 square-degree field from its apex (a fiducial observer 
position) at a corner of the simulation box directed along 
a diagonal of the simulation box up to 2; = 4.9 (Hamana 
et al. 2001b)^. From this light-cone output, we generate 

http://www.mpa- 
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Table 1 

Parameters adopted in the simulation models. 
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VIRGO-SCDM 


1.0 


0.0 


0.0 


0.5 


0.51 


256-^ 


50 


2.27 X 10^" 


0.94 


0.03 


VIRGO-ACDM 


0.3 


0.7 


0.0 


0.7 


0.9 


256^ 


30 


6.86 X 10^° 


0.94 


0.03 


HUBBLE-ACDM 


0.3 


0.7 


0.04 


0.7 


0.9 


109 


35 


2.22 X 10^2 


3.0 


0.1 



three slice subsamples with depth of 200/i~^Mpc centered 
a,t z = 0.5, 1.5 and 2.5. 

We compute the two-point correlation function of dark 
matter particles using the pair-count estimator proposed 
by Landy & Szalay (1993): 

_ DD{r)-2DRir)+RR{r) 

RR{r) • ^^"^^ 



In so doing, we distribute the same number of random 
particles as simulation particles for each model. 




r [h~'Mpc] 

Fig. 3.— Same as Fig. 2 but for VIRGO-ACDM simulation. 




Fig. 4. — Same as Fig. 2 but for Hubble volume ACDM simulation. 



3.2. Model predictions 

We compute model predictions on the basis of the dark 
halo approach which take into account the limitations due 
to the lack of the initial fluctuation power and the finite 
mass resolution by introducing the two parameters fc^ax 
and Mmin as described in §2.3. For definiteness we set 
fcmax as ^Nyquist = tt/ Apart for cach simulation. The min- 
imum mass scale whose density profile can be resolved in 
simulations is determined by the particle mass mpart, but 
the precise value is somewhat uncertain. In the analysis 
below we use M^nin = Afripart with /U = 4, which, we find, 
reasonably fits the A''-body data almost independently of 
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either the redshift or simulation parameters. 

3.3. Evolution of the dark matter clustering at r < Apart 

Figures 2 to 4 present the comparison between the sim- 
ulation data and our model predictions for Virgo SCDM, 
Virgo ACDM and Hubble ACDM data at different red- 
shifts. The filled circles show the two-point correlation 
functions measured directly from each simulation output. 
The mean separation lengths Apart of simulation particles 
are indicated by arrows. Dashed lines show the nonlinear 
prediction by Peacock & Dodds (1996). Predictions by 
the dark halo approach, to which corrections for the lack 
of small-scale perturbation power and the finite particle 
mass are made (see §2), are shown by dotted, long-dashed 
and solid lines for the halo-halo term, Poisson term and 
their sum, respectively. It is important to keep in mind 
that the fitting fmiction by Peacock & Dodds (1996) may 
have 10 to 20 percent uncertainty depending on the cos- 
mological model and scale (Peacock & Dodds 1996), thus 
a small disagreement found at very small scales in lower 
redshift cases may be partly due to a failure of the fitting 
formula (Jenkins et al. 1998). 

At r > Apart) the halo- halo term almost perfectly repro- 
duces the simulation results. This reflects the fact that the 
evolution of large-scale structure is insensitive on average 
to the details of the smaller-scale clustering; incidentally 
this is why we can safely apply linear theory on large scales 
without specifying the small-scale inhomogeneity. 

The evolution of clustering at r < Apart is more inter- 
esting. On these scales, one cannot assign the proper am- 
plitude of fluctuation spectra predicted from linear theory 
at ^Jinit. Nevertheless at z = (3 — 5), the almost right 
amount of power has been already generated by the non- 
linear mode-co\ipling from the larger scales. This suggests 
that the rate of the nonlinear power transfer is sufficiently 
rapid. As a result, at those epochs the small-scale power 
in simulations is primarily limited by the finite mass reso- 
lution rather than the lack of the initial fluctuation power 
on the corresponding scales. 

3.4. When can we trust the small-scale dark matter 

clustering in simulations ? 

It is remarkable that the two-point correlation function 
in simulations on scales below the mean particle separation 
docs not serioiisly sniffer from the lack of the initial fluc- 
tuation power, at least at later epochs. Nevertheless the 
resulting dark matter clustering is still below the nonlinear 
prediction due to the finite mass resolution which cannot 
be accounted within the 10 to 20 percent uncertainty in the 
nonlinear fitting formula (Peacock & Dodds 1996; Jenk- 
ins et al. 1998). In fact, this systematic introduces error 
in the computed correlation function for simulations with 
lower mass resolution and/or at higher redshifts. Then 
one might naturally ask when the small-scale dark matter 
clustering in simulations becomes reliable. We propose an 
empirical criterion as follows. 

We have already shown that the small-scale dark matter 
clustering is primarily limited by the finite mass resolution. 
Also the Poisson term would be dominated by halos which 
become nonlinear at that epoch. The characteristic mass 
for gravitationally collapsed objects, Mnl(-z), is approxi- 
mately given by the condition a{Mf^j,,z) = 1. Therefore 



this mass scale should be much larger than the simulation 
particle mass mpart so as to reproduce the proper ampli- 
tude of the Poisson term. 

Since the small-scale dark matter clustering mainly 
comes from the Poisson term (§2), the absence of halos 
with a mass comparable to and less than the particle mass 
in a A''-body simulation leads to a lack of the clustering 
amplitude at small scales. Therefore, in order to reproduce 
the proper amplitude of the Poisson term at a certain time, 
halos with at least a characteristic mass for gravitation- 
ally collapsed objects at that time must be resolved. Such 
characteristic mass, Mnl(-2), is approximately given by the 
condition (t(Mnl, z) = I, and should be much larger than 
the simulation particle mass mpart so as to be resolved in 
a A^-body simulation. This requirement can be translated 
to the criterion for the critical redshift .Zcrit when simu- 
lations reasonably resolve the dark matter clustering on 
scales below the mean separation: 

MNL(-2crit) = r^halo"^part• (14) 

In the above we introduce a fudge factor nhaio so as to 
represent the number of particles that is required to rea- 
sonably resolve a single halo of the nonlinear mass. Usually 
nhaio is typically supposed to be around 10, but, as far as 
statistical measures such as the dark matter two-point cor- 
relation functions arc concerned, it could be much smaller; 
we used nhaio ~ /x = 4 in plotting Figures 2 to 4, for in- 
stance. 



1 1 1 1 1 1 1 1 1 1 1 1 r 




redshift 

Fig. 5. — The nonlinear mass scale Mnl defined by <jR(Mt^j^) = 1 
versus redshift. The solid and dashed lines are for ACDM and 
SCDM model, respectively. The dotted horizontal lines shows the 
particle masses of AT-body simulations used in §3. The dot-dashed 
line shows the fitting function of MNL-^crit relation for ACDM 
model, equation (15). 

In Figure 5, we plot Mnl(2) against redshift for the 
SCDM (dashed lines) and ACDM (solid lines) parame- 
ters of the present simulations. The particle masses for 
the three simulation models are marked in dotted lines. 
The dot-dashed line shows the approximate relation for 
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the ACDM model that fits the relation for the ACDM 
model within 10 percent accuracy for < 2; < 6; 



= 0.72 



logi 



2 X W^^h-'^Mp 



1.25 



1. 



(15) 



^halo^part 

As far as the dark matter correlations on scales below the 
mean particle separation are concerned, one should not 
trust the simulation data a,t z > Zait- It is important to 
note that this criterion does not mean all results from such 
A'^-body simulations are unreliable, but should be applied 
only to the dark matter correlations at small scales. 

4. SUMMARY AND DISCUSSION 

We have examined the reliability of dark matter corre- 
lation in cosmological high-resolution AT-body simulations 
on scales below the mean particle separation. The parti- 
cle discreteness eflFect imposes two fundamental limitations 
on those scales; the lack of initial fluctuation power and 
finite mass resolution. We applied dark halo approach to 
discuss separately how these two limitations affect the sim- 
ulation results at later epochs. We found that nonlinear 
power transfer below the mean particle separation (up to 
the force resolution) is sufficiently rapid that the reliabil- 
ity of the simulations is primarily determined by the mass 
of particles. By a detailed comparison with three major 
cosmological simulations, we conclude that a,t z > Zait 
the small scale dark matter correlations in high-resolution 
A^-body simulations are not reliable. Of course it is impos- 
sible to prove that the opposite is true at z < 2;crit • Never- 
theless our comparison indicates also that the small scale 
dark matter clustering in a high-resolution A'^-body simu- 
lation is reliable down to the gravitational force resolution 
length for z < Zcnt at least as far as the two-point corre- 
lation functions are concerned. It must be noted that this 
criterion does not mean that all results from such AT-body 
simulations arc unreliable. The limitation of the reliabil- 
ity of a measurement from A^-body simulations depends 
on the quantity one wants to examine. 

The criterion should be applied only to measurements of 
the dark matter clustering such as its two-point correlation 
function and power spectrum. Splinter et al. (1998) con- 
cluded that phase sensitive statistics exhibit strongly the 
effect of discreteness. Strictly speaking our conclusion is 
not inconsistent with their finding, because the two-point 
correlation is insensitive to the phase information and the 
difference would largely come from the higher-order corre- 
lations that we do not attempt to compare. In principle 
our model could be generalized to higher-order correla- 
tions, but there is no reliable empirical model for those 
statistics that we can test against, to the quantitative de- 
gree of the Peacock-Dodds prediction for the two-point 
correlation function. 



We believe that the present methodology to test the 
reliability of the dark matter clustering in a numerical 
simulation is original and powerful, and at least comple- 
mentary to a more traditional way of searching for the 
convergence among different numerical codes with differ- 
ent resolutions. On the other hand, we have to note that 
these approaches are not completely independent; the ma- 
jor ingredients in our model "predictions" including the 
nonlinear mass power spectrum (Peacock & Dodds 1996), 
the mass fiuiction of halos (Shcth & Tormcn 1999), and 
the density profile of halos (Navarro, Frenk & White 1996) 
are actually calibrated by high-resolution simulations. We 
have assumed that those calibrations are performed by per- 
fect simulations, applied the resulting predictions to the 
current problem, and found that our model remarkably 
reproduces what we measure from simulation data. The 
present methodology may still suffer from some calibra- 
tion imccrtainties. Part of this difficulty is already illus- 
trated in the small but systematic deviation of our pre- 
dictions and simulation data on scales one-order-of mag- 
nitude smaller than the mean separation. Most likely this 
is due to the gravitational softening effect in simulation 
data. Nevertheless it is also true that the model predic- 
tion in those scales may be sensitive to the choice of the 
inner slope of halo profile; we adopted the NFW value, 1.0, 
throughout the paper, but it could be 1.5 as many recent 
simulations indicate (Fukushige & Makino 1997; Moore et 
al. 1998). Therefore wc do not claim that the validity of 
the Peacock-Dodds prediction and our halo approach has 
been tested below those scales. Our conclusion, however, 
is valid on scales above ^0.1 times the mean separation 
lengths where those details are not critical. 

Exactly for these reasons, we hope to revisit and re- 
examine the reliability of the model predictions calibrated 
by existing simulations in future. 
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